(* ======================================================================== *)
(*             Basic Multivariate Polynomial Algebra over Q[\vec{x}]        *)
(*                                                                          *)
(* by G.O.Passmore, Cambridge Computer Laboratory and LFCS, Edinburgh, 2011 *)
(* Contact: (e) grant.passmore@cl.cam.ac.uk   (w) www.cl.cam.ac.uk/~gp351/. *)
(* ======================================================================== *)

structure Algebra : Algebra =
struct

(* Var ID type *)

type var_id = int;

(* Var-power type *)

type vp = var_id * int;

(* Power-product type *)

type pp = vp list;

(* Monomial type *)

type mono = Rat.rat * pp;

(* Polynomial type *)

type poly = mono list;

(* Given two power-products in canonical form (ordered with their
   var_id's descending), return their product in canonical form. *)

fun pp_mult' (pp : pp, pp' : pp, result : pp) =
    case (pp, pp') of
	([], a) => (List.rev result) @ a
      | (a, []) => (List.rev result) @ a
      | (((vid, pow) :: s), (((vid', pow') :: s'))) =>
	if (vid = vid') then
	    pp_mult' (s, s', (vid, pow + pow') :: result)
	else if (vid > vid') then
	    pp_mult' (s, pp', (vid, pow) :: result)
	else pp_mult' (pp, s', (vid', pow') :: result);

fun pp_mult (pp : pp, pp' : pp) = pp_mult' (pp, pp', []);

(* Multiply monomials *)

fun m_mult (m : mono, m' : mono) =
    case (m, m') of
	((c, pp), (c', pp')) =>
	(Rat.mult (c, c'), pp_mult(pp, pp'));

(* Negate a monomial *)

fun m_neg (m : mono) = 
    case m of
	(c, pp) => (Rat.neg c, pp);

(* Multivariate total degree of a power-product *)

fun pp_deg' (pp : pp, result) =
    case pp of
	[] => result
      | (_, pow) :: s => pp_deg' (s, pow + result);

fun pp_deg (pp : pp) = pp_deg' (pp, 0);

(* Multivariate total degree of a monomial *)

fun m_deg ((_, p) : mono) = pp_deg p;

(* Deg-diff: The first place in which two monomials disagree *)

fun m_deg_diff' (pp : pp, pp' : pp) =
    case (pp, pp') of
	([], []) => NONE
      | ([], (_, p) :: _) => SOME (0, p)
      | ((_, p) :: _, []) => SOME (p, 0)
      | ((v, p) :: s, (v', p') :: s') =>
	if (v > v') then SOME (p, 0) else
	if (v < v') then SOME (0, p') else
	if not(p = p') then SOME (p, p') else
	m_deg_diff' (s, s');
	    
fun m_deg_diff (m : mono, m' : mono) =
    case (m, m') of
	((_, pp), (_, pp')) =>
	if (pp = pp') then NONE else
	m_deg_diff' (pp, pp');

(* Lexicographic ordering on monomials *)

fun m_lt_lex (m : mono, m' : mono) =
    let val (deg_m, deg_m') = (m_deg m, m_deg m') 
    in
	case m_deg_diff (m, m') of
	    SOME (i, j) => i < j
	  | NONE => false
    end;

(* Graded reverse lexicographic ordering on monomials *)

fun m_lt_degrevlex (m : mono, m' : mono) =
    let val (deg_m, deg_m') = (m_deg m, m_deg m') 
    in
	(deg_m < deg_m') orelse
	((deg_m = deg_m') andalso 
	 case m_deg_diff (m, m') of
	     SOME (i, j) => i > j
	   | NONE => false)
    end;

(* ************************************************************ *)
(* Set monomial order here. Choice is currently between Lex and *)
(* DegRevLex.                                                   *)
(*                                                              *)
(* val m_lt = m_lt_degrevlex;                                   *)
   val m_lt = m_lt_lex;
(* ************************************************************ *)

(* Turn a monomial into a polynomial *)

fun poly_of_mono (m : mono) = 
    case m of 
	(c, _) =>
	if c = Rat.zero then [] else [m];

(* Turn a rational into a polynomial *)

fun poly_of_rat (q : Rat.rat) = poly_of_mono (q, []);

(* Is a monomial's coefficent 0? *)

fun m_zero (m : mono) =
    case m of (c, _) => Rat.eq (c, Rat.zero);

(* Zero polynomial *)

val p_zero = poly_of_rat Rat.zero;

(* One polynomial *)

val p_one = poly_of_rat Rat.one;

(* Is a polynomial a zero polynomial? *)

fun p_isZero p =
    case p of [m] => m_zero m
	    | _ => false;

(* Given two polynomials in canonical form (ordered in descending
   order w.r.t. m_lt), return a new polynomial which is their sum also
   expressed in canonical form. *)

fun p_sum' (p : poly, p' : poly, result : poly) =
    case (p, p') of
	([], s) => if not(p_isZero(s)) then (List.rev result) @ s : poly
		   else (List.rev result)
      | (s, []) => if not(p_isZero(s)) then (List.rev result) @ s
		   else (List.rev result)
      | (m :: s, m' :: s') =>
	case (m, m') of
	    ((c, pp), (c', pp')) =>
	    if (pp = pp') then
		let val d = Rat.add(c, c') in
		    if Rat.eq(d, Rat.zero) then
			p_sum' (s, s', result)
		    else p_sum' (s, s', (Rat.add(c, c'), pp) :: result)
		end
	    else if m_lt(m', m) then
		p_sum' (s, p', if not(m_zero(m)) then m :: result else result) 
	    else p_sum' (p, s', if not(m_zero(m')) then m' :: result else result)
	    
fun p_sum (p : poly, p' : poly) = p_sum' (p, p', [] : poly);

(* Negate a polynomial in canonical form *)

fun p_neg (p : poly) = List.map m_neg p;

(* Subtract two polynomials in canonical form *)

fun p_sub (p : poly, p' : poly) =
    p_sum (p, p_neg p') : poly;

(* Multiply a monomial and a polynomial, both in canonical form *)

fun mp_mult' (m : mono, p : poly, result : poly) =
    case p of
	[] => result
      | (m' :: s') =>
	let val x = m_mult(m, m') in
	    mp_mult' (m, s', p_sum(poly_of_mono(x), result))
	end;

fun mp_mult (m : mono, p : poly) = mp_mult' (m, p, []) : poly;

(* Multiply two polynomials, both in canonical form *)

fun p_mult' (p : poly, p' : poly, result : poly) =
    case (p, p') of
	([], _) => result
      | (m :: s, p') =>
	p_mult'(s, p', p_sum(result, mp_mult(m, p')));

fun p_mult (p : poly, p' : poly) =
    if List.length(p) <= List.length(p') then
	p_mult' (p, p', []) else p_mult' (p', p, []) : poly;

(* Make a polynomial monic (used when it is implicitly = 0) *)

fun p_monic (p : poly) =
    let val (c, _) = List.hd p
	val c' = ((Rat.inv c), []) : mono
    in mp_mult (c', p) end;

(* Compute the LCM of two power-products in canonical form *)

fun pp_lcm' (pp : pp, pp' : pp, result : pp) : pp = 
    case (pp, pp') of
	([], pp') => (List.rev result) @ pp'
      | (pp, []) => (List.rev result) @ pp
      | ((v, pow) :: r, (v', pow') :: r') =>
	if (v > v') then
	    pp_lcm' (r, pp', (v, pow) :: result) 
	else if (v < v') then
	    pp_lcm' (pp, r', (v', pow') :: result)
	else if (pow > pow') then
	    pp_lcm' (r, r', (v, pow) :: result)
	else pp_lcm' (r, r', (v', pow') :: result);

fun pp_lcm (pp : pp, pp' : pp) = pp_lcm' (pp, pp', []) : pp;

(* Divide a pp by another, breaking if division would be unclean. 
   Again, we utilise the fact that the pps are in canonical form. *)

fun pp_div' (pp : pp, pp' : pp, result : pp) : pp =
    case (pp, pp') of
	([], []) => (List.rev result)
      | ([], _) => raise Useful.Error "Unclean power-product division"
      | (pp, []) => (List.rev result) @ pp
      |  ((v, pow)::r, (v', pow')::r') =>
	 if v' > v then raise Useful.Error "Unclean power-product division"
	 else if v > v' then pp_div' (r, pp', (v, pow) :: result)
	 else if pow > pow' then pp_div' (r, r', (v, (pow - pow')) :: result)
	 else if pow = pow' then pp_div' (r, r', result)
	 else raise Useful.Error "Unclean power-product division";

fun pp_div (pp : pp, pp' : pp) : pp = pp_div' (pp, pp', [] : pp) : pp;

(* Divide one monomial by another, breaking if the division
   would be unclean. *)

fun m_div (m : mono, m' : mono) =
    case (m, m') of
	((c, pp), (c', pp')) =>
	(Rat.mult (c, Rat.inv(c')),
	 pp_div (pp, pp'));   

(* Does one power-product, pp, divide another one, pp'? *)

fun pp_divides (pp : pp, pp' : pp) =
    case (pp, pp') of
	([], _) => true
      | (_, []) => false
      | ((v, pow)::r, (v', pow')::r') =>
	if (v > v') then false else
	if (v = v') then 
	    (pow <= pow' andalso pp_divides (r, r')) 
	else pp_divides (pp, r');  

(* Given a polynomial, return its head monomial. *)

fun p_hm (p : poly) : mono = 
    case p of 
	[] => (Rat.zero, []) 
      | (c, pp) :: _ => (c, pp);

(* Given a polynomial, return its head power-product. *)

fun p_hpp (p : poly) : pp = case p of [] => [] | (_, pp) :: _ => pp;

(* Are two monomials equal? *)

fun m_eq (m, m') =
    case (m, m') of
	((c, pp), (c', pp')) =>
	Rat.eq (c, c') andalso pp = pp';

(* Are two polynomials equal? *)

fun p_eq (p, p') =
    List.length p = List.length p'
    andalso ListPair.all m_eq (p, p');

(* Power-product to string *)

fun pp_toString (pp : pp) =
    case pp of
	[] => ""
      | ((v, p) :: s) => "x" ^ Int.toString(v)
			 ^ (if not(p = 1) then ("^" ^ Int.toString(p)) else "")
			 ^ (if null(s) then "" else (" " ^ pp_toString(s)));
	
(* Monomial to string *)

fun m_toString (m : mono) =
    case m of
	(c, pp) => 
	let val pp_null = null(pp) in
	    (if Rat.eq(c, Rat.one) andalso not(pp_null) then "" 
	     else (Rat.toString(c) ^ (if pp_null then "" else " ")))
	    ^ pp_toString(pp)
	end;

(* Polynomial to string *)

fun p_toString (p : poly) =
    case p of
	[] => ""
      | (m :: s) => m_toString(m) ^ 
		    (if null(s) then "" else (" + " ^ p_toString(s)));

end
